Carregando os dados

df <- read_excel("DADOS_BR_SB.xlsx")
# Banco do Brasil SA (BBAS3)
bbas <- na.omit(xts(as.numeric(df$BBAS3_SB), order.by=as.Date(df$Date)))
bbas <- na.omit(100*diff(log(bbas)))
# Bovespa Index (IBOV)
ibov <- na.omit(xts(as.numeric(df$IBOVC_SB), order.by=as.Date(df$Date)))
ibov <- na.omit(100*diff(log(ibov)))
# Selic 252
selic <- na.omit(xts(as.numeric(df$Selic252_SB), order.by=as.Date(df$Date)))
selic <- na.omit(100*diff(log(selic)))

Questão 1

1.

O modelo CAPM que o excesso de retornos do BBAS3 e do mercado (IBOV) são: Estocástico estacionário, serialmente não correlacionado e com distribuição normal.

Vamos analisar essas hipóteses uma a uma. Primeiro vamos gerar a série dos nossos excessos de retornos.

# Excesso de retorno BBAS3
exrcp_bbas <- bbas - selic
plot.xts(exrcp_bbas)

# Excesso de retorno do mercado (IBOV)
exrcp_ibov <- ibov - selic
plot.xts(exrcp_ibov)

Vamos começar pelo Excesso de retorno BBAS3, analisando a estacionaridade.

# ACF and PACF
par(mfrow=c(1,2))
Acf(exrcp_bbas, lag.max=45, main = "IBOV")
Pacf(exrcp_bbas, lag.max=45, main = "IBOV")

# Augmented Dickey-Fuller
adf.test(exrcp_bbas)

## Warning in adf.test(exrcp_bbas): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  exrcp_bbas
## Dickey-Fuller = -17.439, Lag order = 18, p-value = 0.01
## alternative hypothesis: stationary

ndiffs(exrcp_bbas,test="adf")

## [1] 0

A série do excesso de retornos do BBAS3 parece não ser estacionária dada nossa análise da FAC e FACP. O teste ADF no entanto indica que a série é estacionária e nenhuma transfomarção é precisa para induzir estácionariedade.

Pela FAC também temos que nossa série apresenta correlação serial, violando a hipótese de que os excessos de retornos do BBAS3 sejam serialmente não correlacionado.

Agora vamos analisar a normalidade.

# Histograma
m.bbas = mean(exrcp_bbas, na.rm=TRUE)
std.bbas = sqrt(var(exrcp_bbas, na.rm=TRUE))
par(mfrow=c(1,1))
hist(exrcp_bbas, breaks=100, prob=T, 
     xlab="BBAS3", ylim=c(0, .25), xlim=c(-15,15),
     main="Excesso de retorno BBAS3")
lines(density(exrcp_bbas), col = "red")   
curve(dnorm(x, mean=m.bbas, sd=std.bbas), 
      col="darkblue", lwd=2, add=TRUE, yaxt="n")

# Jarque–Bera
jarque.bera.test(exrcp_bbas)

## 
##  Jarque Bera Test
## 
## data:  exrcp_bbas
## X-squared = 34453, df = 2, p-value < 0.00000000000000022

Pelo teste de Jarque–Bera vemos que a nossa série não é normal e logo a hipótese de normalidade não é valida para o excesso de retornos do BBAS3.

Agora vamos analisar o excesso de retorno do mercado (IBOV), começando pela estacionaridade.

# ACF and PACF
par(mfrow=c(1,2))
Acf(exrcp_ibov, lag.max=45, main = "IBOV")
Pacf(exrcp_ibov, lag.max=45, main = "IBOV")

# Augmented Dickey-Fuller
adf.test(exrcp_ibov)

## Warning in adf.test(exrcp_ibov): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  exrcp_ibov
## Dickey-Fuller = -15.894, Lag order = 18, p-value = 0.01
## alternative hypothesis: stationary

ndiffs(exrcp_ibov,test="adf")

## [1] 0

Assim como no excesso de retornos do BBAS3, a série do excesso de retornos do IBOV parece não ser estacionária dada nossa análise da FAC e FACP. O teste ADF no entanto indica que a série é estacionária e nenhuma transfomarção é precisa para induzir estácionariedade.

Pela FAC também temos que nossa série apresenta correlação serial, violando a hipótese de que os excessos de retornos do IBOV sejam serialmente não correlacionado.

Analisando a normalidade.

# Histograma
m.ibov = mean(exrcp_ibov, na.rm=TRUE)
std.ibov = sqrt(var(exrcp_ibov, na.rm=TRUE))
par(mfrow=c(1,1))
hist(exrcp_ibov, breaks=100, prob=T, 
     xlab="IBOV", ylim=c(0, .35), xlim=c(-10,10),
     main="Excesso de retorno IBOV")
lines(density(exrcp_ibov), col = "red")   
curve(dnorm(x, mean=m.ibov, sd=std.ibov), 
      col="darkblue", lwd=2, add=TRUE, yaxt="n")

# Jarque–Bera
jarque.bera.test(exrcp_ibov)

## 
##  Jarque Bera Test
## 
## data:  exrcp_ibov
## X-squared = 273736, df = 2, p-value < 0.00000000000000022

Novamente pelo teste de Jarque–Bera vemos que a nossa série não é normal e logo a hipótese de normalidade não é valida para o excesso de retornos do IBOV.

Com base em nossa análise temos que essas hipóteses não são adequadas pois as séries não são normais nem serialmente não correlacionado.

2. Estimação por MQO

capm <- lm(exrcp_bbas ~ exrcp_ibov)
summary(capm)

## 
## Call:
## lm(formula = exrcp_bbas ~ exrcp_ibov)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -19.9965  -1.0230  -0.0631   0.9681  11.2875 
## 
## Coefficients:
##             Estimate Std. Error t value            Pr(>|t|)    
## (Intercept)  0.03388    0.02415   1.403               0.161    
## exrcp_ibov   1.05728    0.01131  93.511 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.851 on 5877 degrees of freedom
## Multiple R-squared:  0.5981, Adjusted R-squared:  0.598 
## F-statistic:  8744 on 1 and 5877 DF,  p-value: < 0.00000000000000022

Por violar duas das hipóteses acima, os parametros estimados por MQO não são consistentes. Temos \(\alpha\) é não significativo, o que implica que o retorno do BBAS3 correspondeu ao benchmark (IBOV) no período.

3. Estimação por MQR

## Estimando Recursivamente
df_recursivo <- cbind(exrcp_bbas,exrcp_ibov)
colnames(df_recursivo) <- c("exrcp_bbas", "exrcp_ibov")
eq1r <- lapply(seq(10, length(exrcp_bbas)), function(x) lm(exrcp_bbas ~ exrcp_ibov, data = df_recursivo[1:x, ]))
# Obtendo todos os coeficients
all_intercepts <- unlist(sapply(1:length(eq1r),function(j) eq1r[[j]]$coefficients[1]))
all_slopes <- unlist(sapply(1:length(eq1r),function(j) eq1r[[j]]$coefficients[2]))
# sd
sd_intercepts <- unlist(sapply(1:length(eq1r),function(j) summary(eq1r[[j]])$coefficients[3]))
sd_slopes<-unlist(sapply(1:length(eq1r),function(j) summary(eq1r[[j]])$coefficients[4]))
# Plotando coeficientes
par(mfrow = c(1,1))
## Alfa
plot.zoo(cbind(all_intercepts, 
               all_intercepts+ sd_intercepts*1.96, 
               all_intercepts-sd_intercepts*1.96), 
         plot.type = "single", 
         main = "Estimação Recursiva - Alfa", 
         ylab = "", xlab = "",
         col = c(1,2,2))

## Beta
plot.zoo(cbind(all_slopes, 
               all_slopes + sd_intercepts*1.96, 
               all_slopes - sd_intercepts*1.96), 
         plot.type = "single", 
         main = "Estimação Recursiva - Beta", 
         ylab = "", xlab = "",
         col = c(1,2,2))

O parametro \(\alpha\) tem valor proximo ao estimado por MQO e a janela estimada por MQR está dentro da estimativa por MQO. O parametro \(\beta\) pode estar superestimado por MQO ja que a estimativa recursiva cresce ligeralmente ao longo do tempo e é menor do que 1.0 no início da amostra.

4.

Estimando por o modelo

\[ y_t = \alpha_t +\beta_t \cdot x_t \] \[ \beta_t = \beta_{t-1} + \eta_t \] \[ \alpha_t = \bar{y_t} - \beta_t \cdot \bar{x_t} \]

model.tvLM <- tvLM(exrcp_bbas~exrcp_ibov)

## Calculating regression bandwidth... bw =  0.06121848

summary(model.tvLM)

## 
## Call: 
## tvLM(formula = exrcp_bbas ~ exrcp_ibov)
## 
## Class:  tvlm 
## 
## Summary of time-varying estimated coefficients: 
## ================================================ 
##         (Intercept) exrcp_ibov
## Min.      -0.147737     0.5118
## 1st Qu.   -0.007129     0.9527
## Median     0.032933     1.0201
## Mean       0.028401     1.0843
## 3rd Qu.    0.067817     1.1205
## Max.       0.195757     1.6659
## 
## Bandwidth:  0.0612
## Pseudo R-squared:  0.622

# 90% IC do Beta
model.tvLM.90 <- confint(model.tvLM, level = 0.95, runs = 50)
plot(model.tvLM.90)

Vermelho = MQO

Preto = MQR

Azul = estimação variável

# Alfas
a <- seq(1,5879) |> as_tibble()
a$variavel <- model.tvLM$coefficients[,1]
a$recursivo <- c(all_intercepts,rep(NA,9))

ggplot(a,aes(x=value,y=recursivo)) +
  geom_line() +
  geom_line(aes(y=variavel), color= "blue") +
  geom_hline(yintercept=0.03388,color = "red") +
  labs(x="",y="",title = "Estimação do Alfa")

## Warning: Removed 9 row(s) containing missing values (geom_path).

Analisando os Alfas, vemos que o Alfa estimado por MQO está subestimado no ínicio da amostra e sobre estimado no final da amostra em relação a estimação variável. O Alfa estimado por MQR apresenta a mesma trajetoria de crescimento no inicio da amostra, porem ele fica sobre estimado no meio ao final da amosta em relação a estimação variável.

# Betas
a <- seq(1,5879) |> as_tibble()
a$variavel <- model.tvLM$coefficients[,2]
a$recursivo <- c(all_slopes,rep(NA,9))

ggplot(a,aes(x=value,y=recursivo)) +
  geom_line() +
  geom_line(aes(y=variavel), color= "blue") +
  geom_hline(yintercept=1.05728, color = "red") +
  labs(x="",y="",title = "Estimação do Beta")

## Warning: Removed 9 row(s) containing missing values (geom_path).

Analisando os Betas vemos que a estimação por MQO sempre apresenta um resultado subótimo em relação aos outros dois métodos. O interessante é que o o MQR parece sub estimar o Beta em relação ao modelo com estimação variável na maior parte da amostra.

Questão 2

# Dados Mensais
df$Date <- as.Date(df$Date , format="%Y-%m-%d" )
df$Ano <- format(as.Date(df$Date, format="%Y/%m/%d"),"%m")
Date <- c(NA,NA)
for (i in 1:(nrow(df)-1)){
  if (df$Ano[i] != df$Ano[i+1]){
    Date[i] <- format(as.Date(df$Date[i], format="%Y-%m-%d"),"%Y-%m-%d")
}}

Date <- Date[is.na(Date) == FALSE]
df_m <- data.frame(Date)
df_m$Date <-as.Date(df_m$Date , format="%Y-%m-%d" )

# Preço Mensal
df_m$bbas_m <- NA
df_m$ibov_m <- NA
for ( i in 1:nrow(df_m)){
  df_m$bbas_m[i] <- df$BBAS3_SB[which(df$Date == df_m$Date[i])]
  df_m$ibov_m[i] <- df$IBOVC_SB[which(df$Date == df_m$Date[i])]
}

# Retorno Mensal
df_m$lbbas_m <- NA
df_m$libov_m <- NA
for ( i in 2:nrow(df_m)){
  temp_bbas = na.omit(100*diff(log(df_m$bbas_m)))
temp_libov = na.omit(100*diff(log(df_m$ibov_m)))
  df_m$libov_m[i] <- temp_libov[i]
  df_m$lbbas_m[i] <- temp_bbas[i]
}

# Retorno Mensal como Time Series
libov_m = ts(df_m$libov_m[2:269],start=c(2000,1),frequency=12)
lbbas_m = ts(df_m$lbbas_m[2:269],start=c(2000,1),frequency=12)
# Preço Mensal como Time Series
ibov_m = ts(df_m$ibov_m[1:270],start=c(2000,1),frequency=12)
bbas_m = ts(df_m$bbas_m[1:270],start=c(2000,1),frequency=12)

Sazonalidade Ibov

# Plotando a série do log retorno do ibov assim como os testes ACF e PACF
plot(libov_m)

Acf(libov_m, lag.max=60, main ='')

Pacf(libov_m, lag.max=60, main ='')

A análise das FAC e FACP da primeira diferença do Log IBOV não nos dá fortes indicios de presença de sazonalidade na série. Na FAC podemos observar alguns pontos saindo da banda. Na FACP observamos apenas um ponto fora da banda, entre o lag 24 e o lag 36. No entanto, temos que lembrar que essa banda é menor do que a banda correta. Então vamos utilizar outros métodos exploratórios para avaliar uma possível sazonalidade em nossa série.

Vamos utilizar a função decompose(), que através de médias móveis nos dá os componentes sazonal, trend e irregular do IBOV.

# Decompondo a série do Ibov
plot(decompose(libov_m,type = "additive"))

plot(decompose(libov_m,type = "multiplicative"))

O modelo “additive” usado é: \[Y_t = T_t + S_t + e_t \]

O modelo “multiplicative” usado é:

\[Y_t = T_t S_t e_t \]

Em ambos os tipos de decomposição, “additive” e “multiplicative”, quando analisamos o gráfico do componente sazonal da série, a série do Ibov aparenta ter um comportamento sazonal determinístico.

Seguimos para uma segunda análise do componente sazonal da nossa série usando a função monthplot(). Essa função extrai subséries mensais de uma série temporal e plotam com suas médias.

# Subséries do Ibov
monthplot(libov_m, col.base = 2, lty.base = 2)

A análise das subséries mensais nos indica a presença de um componente sazonal estocástico. Vemos que as médias não são muito diferentes para os meses, no entanto temos que a variância dentro dos meses não é constante.

A seguir usamos a função seas() para realizar um ajuste sazonal X-13ARIMA-SEATS em nossa série.

ajuste_ibov <- seas(libov_m)

## Model used in SEATS is different: (0 0 0)

summary(ajuste_ibov)

## 
## Call:
## seas(x = libov_m)
## 
## Coefficients:
##                    Estimate Std. Error z value    Pr(>|z|)    
## AO2020.Jan        -35.59708    6.79605  -5.238 0.000000162 ***
## MA-Nonseasonal-01  -0.15877    0.06074  -2.614     0.00895 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## SEATS adj.  ARIMA: (0 0 1)  Obs.: 268  Transform: none
## AICc:  1801, BIC:  1811  QS (no seasonality in final):0.9008  
## Box-Ljung (no autocorr.): 19.02   Shapiro (normality): 0.9902 .
## Messages generated by X-13:
## Warnings:
## - Automatic transformation selection cannot be done on a series
##   with zero or negative values.
## - Series should not be a candidate for seasonal adjustment
##   because the spectrum of the prior adjusted series (Table B1)
##   has no visually significant seasonal peaks.
## 
## Notes:
## - Model used for SEATS decomposition is different from the model
##   estimated in the regARIMA modeling module of X-13ARIMA-SEATS.

plot(ajuste_ibov)

A transformação automática não pode ser realizada em nossa série pois ela não é estritamente positiva (Vemos que a série ajustada pelo seas() é igual a série original).

# QS statistic
qs(seas(libov_m))

## Model used in SEATS is different: (0 0 0)

##                   qs   p-val
## qsori        0.90084 0.63736
## qsorievadj   2.12199 0.34611
## qsrsd        2.30881 0.31525
## qssadj       0.90084 0.63736
## qssadjevadj  2.12199 0.34611
## qsirr        1.23580 0.53908
## qsirrevadj   2.19019 0.33451
## qssori       0.59092 0.74419
## qssorievadj  5.50828 0.06366
## qssrsd       4.08100 0.12996
## qsssadj      0.59092 0.74419
## qsssadjevadj 5.50828 0.06366
## qssirr       1.18405 0.55321
## qssirrevadj  4.12419 0.12719

O teste de sazonalidade da função seas() é dado pela estatística QS e a hipótese nula do teste é não-sazonalidade. Nós não rejeitamos a hipótese nula com 95% de nivél de confiânça.

Sazonalidade BBAS3

# Plotando a série do log retorno do BBAS3 assim como os testes ACF e PACF
plot(lbbas_m)

Acf(lbbas_m, lag.max=60, main ='')

Pacf(lbbas_m, lag.max=60, main ='')

A análise das FAC e FACP da primeira diferença do Log BBAS3 não indica presença de sazonalidade na série. Na FAC temos poucos pontos saindo da banda (novamente, a banda correta é maior) e na FACP não temos nenhum ponto saindo da banda depois do 12th lag. Assim como anteriormente, isso indica que a série provavelmente não tem um comportamento média-móvel sazonal. Vamos seguir com outros métodos de análise.

# Decompondo a série do BBAS3
plot(decompose(lbbas_m,type = "additive"))

plot(decompose(lbbas_m,type = "multiplicative"))

Analisando os gráficos do componente sazonal da função decompose() temos que aparentemente nossa série possuí um comportamento sazonal determinístico. Podemos observar que a variância do componente sazonal é constante em todo o período análizado.

# Subséries do BBAS3
monthplot(lbbas_m, col.base = 2, lty.base = 2)

A análise do gráfico da função monthplot() nos sugere a presença de um componente sazonal estocástico. Temos choques de volatilidade em alguns meses como Janeiro, Março, Junho e Julho. Ou seja, nosso componente sazonal não parece ter variância constante. Mas, as médias são muito parecidas em todos os meses.